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Abstract 

We present a method for estimating the density of states of a classical 
statistical model. The algorithm successfully combines the Wang-Landau 
flat histogram method with the N-fold way in order to improve efficiency 
of the original single spin flip version. We test our implementation of the 
Wang-Landau method with the two-dimensional nearest neighbor Ising 
model for which we determine the tunneling time and the density of states 
on lattices with sizes up to 50 x 50. Furthermore, we show that our new al- 
gorithm performs correctly at right edges of an energy interval over which 
the density of states is computed. This removes a disadvantage of the 
original single spin flip Wang-Landau method where results showed sys- 
tematically higher errors in the density of states at right boundaries. In 
order to demonstrate the improvements made, we compare our data with 
the detailed numerical tests presented in a study by Wang and Swendsen 
where the original Wang-Landau method was tested against various other 
methods, especially the transition matrix Monte Carlo method (TMMC). 
Finally, we apply our method to a thin Ising film of size 32 x 32 x 6 with 
antisymmetric surface fields. With the density of states obtained from the 
simulations we calculate canonical averages related to the energy such as 
internal energy, Gibbs free energy and entropy, but we also sample mi- 
crocanonical averages during simulations in order to determine canonical 
averages of the susceptibility, the order parameter and its fourth order 
cumulant. We compare our results with simulational data obtained from 
a conventional Monte Carlo algorithm. 



1 Introduction 

While Monte Carlo methods in statistical thermodynamics already find broad 
application , the standard approach using the Metropolis algorithm Q has 
the important disadvantage that the entropy of the simulated model system 
is not an output of the calculation. Furthermore in systems with a rugged 
landscape of the (coarse-grained) free energy the convergence of the method in 
practice often is problematic - the system may get "stuck" in one valley of this 
free energy landscape for a long time. 

In order to overcome these difficulties, many interesting approaches (e.g., um- 
brella sampling || , multicanonical Monte Carlo || , expanded ensemble methods 
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(7[||) have been proposed in the literature, and some of these techniques em- 
phasize the idea of directly sampling the energy density of states (e.g., 
Such approaches are clearly very promising, and in particular the method of 
Wang and Landau Jl2j has the additional merit that it is straightforward to im- 
plement. However, it is a difficult matter to judge the efficiency of the various 
methods, and estimate in beforehand the computational effort that is necessary 
to reach a desired level of accuracy for a chosen linear dimension L of the (lat- 
tice) model under study. In particular, a recent comparative study of Wang 
and Swendsen [fl4|| cautioned against a too optimistic view on these matters, 
and suggested that the Wang-Landau method 12 for medium and large lattice 
sizes L is considerably less efficient than various versions of transition matrix 
Monte Carlo methods [fl4|| , which are more complicated to implement however. 
In the present paper, we reconsider the Wang-Landau algorithm Jl^] and modify 
it by combining it with the N-fold way algorithm of Bortz et al. 15 1. It has been 
known that this latter algorithm is far superior to the single spin flip Metropolis 
algorithm |l| at low temperatures. However, as will be demonstrated here, our 
new algorithm performs significantly better near the right edge of an energy 
interval, since for large L it is necessary to split the total energy range in many 
subintervals (which are slightly overlapping, of course, so that they can be joined 
unambigously) , which then can be sampled in parallel on the processors of a 
multi-processor machine. This improvement of accuracy on the right boundaries 
of energy intervals is crucial for obtaining a very good overall accuracy. Thus, 
we are able to show that this particular criticism does no longer apply with 
respect to this new algorithm and a performance of the algorithm comparable 
in efficiency to the transition matrix Monte Carlo methods is reached. 
In sec. |2|, the flat histogram method of Wang and Landau Jl2| is briefly re- 
viewed, as well as the N-fold way algorithm |l5), and our combination of these 
two methods is presented. Sec. ^ then describes the implementation and re- 
sults for the two-dimensional Ising model, treated also in refs. ]I^-|l4j|, and a 
comparative analysis of errors in the density of states is presented. Sec. [I] then 
gives, as an example of an application to a nontrivial model of current interest 
|l6|-[l8[ results for a three-dimensional Ising model in a L x L x D thin film 
geometry (L — 32, D = 6), where competing boundary fields (hi = — /12) act 
on the two free L x L surfaces. It is shown that the present method produces 
results compatible with those of the conventional single spin flip algorithm [filif 
and reaches a better accuracy with less effort in computing time. Sec. || finally 
summarizes some conclusions. 



2 The flat histogram method of Wang-Landau 
2.1 Single spin flip 

Recently, F. Wang and Landau |l^] proposed a Monte Carlo algorithm for clas- 
sical statistical models which uses a random walk in energy space in order to 
obtain an accurate estimate for the density of states g(E). This method is based 
upon the fact that a flat energy histogram H(E) is produced if the probability 
for the transition to a state of energy E is proportional to l/g(E). 
This observation is utilized in the following way. Initially, g(E) is set equal to 
one for all energies. A spin is then chosen at random and flipped with probabil- 
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ity min(l, g(E)/g(E')) whereby E' is the energy of the system with the chosen 
spin being overturned. The density of states g(E) is not constant during the 
random walk, but is updated according to g(E) — > g{E) ■ f after each spin flip 
trial whether the spin is flipped or not. A histogram H(E) records how often a 
state of energy E is visited. In the beginning of the random walk the modifica- 
tion factor / can be as large as e ~ 2.7182818.[|Each time the energy histogram 
satisfies a certain flatness criterion, / is reduced according to / — > /3 and H(E) 
is set to zero for all energies. The histogram is considered as flat if 

H(E)>e-{H(E)) (f) 

for all E where e is usually between 0.7 and 0.95 and (H(E)) is the average 
histogram. The simulation is ended if / is close enough to one, i.e., smaller than 
some predetermined final modification factor f final- To speed up simulations 
it is possible to perform several random walks on adjacent energy intervals on 
independent processors. Disadvantages of the single spin flip version of this 
algorithm are the small acceptance rates for energy intervals which contain the 
groundstate and low- lying excited states and the relatively large errors of g(E) 
at right edges of energy intervals as reported in (l2) . Since this flat histogram 
method produces only a relative density of states g(E), one has to normalize 
g(E) in order to get the absolute density of states g(E). This can be done 
by using known values of the density of states, for example, the groundstate 
degeneracy or other constraints on g(E). In case of the two-dimensional nearest 
neighbor Ising model 

Ti = — J SiSj with Si — ±1, (2) 

one has 



g(-2JN) = 2, (3) 
g(-E) = g(E), (4) 

Es(£) = 2Ar > ( 5 ) 

E 

where N is the number of spins. If one uses the groundstate to normalize 
g(E) one can calculate canonical averages from the density of states for all 
temperatures which become exact for T — > 0. It is thus highly desirable to 
estimate the density of states for low-lying energies with sufficient accuracy. 
Since we deal with low acceptance rates of the single spin flip approach in these 
energy-ranges one can reduce the simulational effort enormously by using a 
rejection-free update scheme. In jlj] J. S. Wang and Swendsen presented an 
efficient Monte Carlo method using the N-fold way to obtain an estimate for 
the transition matrix from which the density of states can be determined via 
an optimization procedure. They have tested their algorithm with the two- 
dimensional Ising model and compared it to various other methods from which 
the density of states can be calculated. In particular, they found their method 
to perform better than the single spin flip Wang-Landau flat histogram method. 

1 If one chooses to sample S'(-E) = logg(E) the modification factor becomes an increment 
S(E) -> S(E) + log / with log / < loge ~ 0.4342944. 
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Adopting the idea of Wang and Landau in the context of transition matrix they 
improved the efficency of the original algorithm only for small system sizes up 
to L = 8, but for larger system sizes they reported that one has the problem of 
sticking to a Gaussian distribution for the histogram. In the next subsection we 
will discuss how the algorithm of Wang and Landau can be successfully combined 
with the rejection-free N-fold way in order to enhance its performance. 



2.2 N-fold way 

In the N-fold way @] a flip occurs at each step of the algorithm and one 
then calculates the life-time of the resulting state. The usual single spin flip 
dynamics is thus preserved and observables become life-time weighted averages 
over the generated states. We will now describe the method in the context of 
the algorithm of Wang and Landau. Since the density of states g(E) can be 
very large especially for large system sizes we consider S(E) — log g(E) during 
simulations. In the beginning S(E) is set to zero for all E. Initially the system 
may be in the state a with energy E G I = [E m i ni E max ], whereby I denotes 
the energy-range for which one wants to estimate g(E). One then partitions 
all spins into classes according to their energetic local environment, i.e., the 
energy difference AEi a spin flip will cause. For a two-dimensional nearest 
neighbor Ising model each spin belongs to one of only M = 10 classes. The 
total probability P of any spin of class i being overturned is given by 

P(AE i ) = n(a,AE i )p(E^E + AE i ), i = l,...,M, (6) 

with n(a, AEi) being the number of spins of state a which belong to class i and 
p(E — > E + AEi) is given by 

n(F^ F+AF\- f min(l, <?(£)/<?(£ + if E + AE, e I 

p{E^E + AE t )-^ Q . f E + AE .g L (7) 

In order to determine the class from which to flip a spin one calculates the 
numbers 

Q m = ^P{AE l ), m = l,...,M and Q = 0, (8) 

which are the integrated probabilities for a spin flip within the first m classes. 
Hence a class is selected by generating a random number < r < Qm and 
taking class m if Q m ~i < r < Q rn . The spin to be overturned is chosen from 
this class with equal probabilities. Due to the flip, the spin and its interacting 
neighbors will change classes and correspondingly the numbers n(cr, AEi) will 
differ from their predecessors. Finally, one has to determine the average life-time 
r of the resulting state, i.e., one has to find out how many times the move just 
made would be rejected on average in the usual update scheme. The probability 
that the first random number would produce a flip is P = Qm/N. Therefore 
one has for the probability that exactly n random numbers will result in a new 
configuration 

P n =P(l-P) n -\ (9) 
Thus the average life-time becomes 



EnP„^ nP(l - P)"" 1 = * (10) 
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Now, the steps which are actually carried out by the N-fold way version of the 
Wang-Landau method are the following: 

1. choose an initial configuration and set H(E) = 0, S(E) = for all E and 
/o = loge ~ 0.4342944 and also fix ff ina i. 

2. determine (update) the probabilities p(E — > E + AEi) and the Q m 's of 
the (initial) configuration using eqs. (||), (0) and (g). 

3. determine average lifetime r of (initial) configuration via eq. (|l0|). 

4. increment histogram, density of states and update ff. 

H{E) -> H(E) + t 
S{E) -» S(£) + AS(£) 

with: 

AS(E) = ( * loge -0.4342944 

v ; (_ log e if /, • t > log e y ' 

f J /i if ji-T< loge 

\ AS(£)/r if / 4 -r > loge, ^ 

in case of fo/fi+i > 2 we set -> /j/2. 

5. after some fixed sweeps check H(E) and refine /j according to = fj/2 
if H(E) is flat. Stop if histogram is flat for a /j < f final- 

6. determine the Q m 's (the p(i? — > £J + A£Jj)'s are not updated here). 

7. choose and flip spin as described above. 

8. go to 2. 



Equations ( |ll| ) and (|12|) ensure that the increment AS(E) is kept below or equal 
to loge. In the single spin flip case this constraint is always guaranteed by the 
choice of the initial modification factor. Larger upper limits for AS(E) result 
in large statistical errors in g(E) as was reported in Jl2| . Furthermore they can 
lead to Qm = in the context of the N-fold way which will definitely destroy 
the iteration procedure described above. We have observed that adjusting 
according to eqs. ( p"l| ) and (|l2|) only applies to the first stage of iteration. Once 
the histogram is flat for the first time the increment AS(E) is simply /j • r. 
Controling /j in other ways may also work. 



3 Simulations 

3.1 Two-dimensional Ising model 

First we have tested our algorithm with the two-dimensional Ising model, eq. 
(||), on L x L square lattices with L = 8, 16, 32, 50 over the entire energy-range 
E/JN £ [—2,2]. For these system sizes, the density of states g(E) can be 
calculated exactly p0| and we will denote it by g ex (E). As already mentioned 
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the method only provides a relative density of states. To obtain the absolute 
density of states g{E), we first utilize ([|) by taking 



g(E)^(g(E)+g(-E))/2. 



(13) 



In order to normalize g(E) one can now use eq. which will yield the best 
accuracy in the density of states for low and high-lying energy levels or one 
uses eq. (||) which will lead to best accuracy in g(E) for energies for which the 
corresponding density of states has large values. We have to stress here that this 
is the simplest method to obtain an absolute density of states, but it suffices in 
order to demonstrate the improvements made by combining the Wang-Landau 
method with the N-fold way. For further improvement of accuracy concerning 
the density of states, one has to resort to optimization methods. This was 
done succesfully in jl4| , where a least-squares method was considered with eqs. 
(||),(§|) and (||) as constraints. 

To compare our data with the tests provided in [jl4| we consider the relative 
error per energy level for the density of states 



s(E) 



9(E) 

9ex{E) 



(14) 



and its average 

E 

The factor 1 / (N — 2) is due to the fact that we have N — 1 different energy levels 
whereby the groundstate is excluded from averaging because we have imposed 
the exact value for the corresponding density of states, see eq. (|J). For each 
system size we have performed 30 runs, starting always with all spins up. Monte 
Carlo time was measured in units of sweeps, whereby one sweep was taken to 
be N spin flips. The results are given in table |l| and figures || and ||. The 
used parameters e and f final are stated in the corresponding captions. We also 
measured the tunneling-time Tt, i.e., the average number of sweeps which is 
elapsing during a transition from the groundstate to a state of highest energy 
or vice versa. The tunneling time is averaged over the complete run, i.e., the 
various stages with different modification increments. For the two-dimensional 
Ising model we get 

r t si 0.79L 2 58 (16) 

from a fit which is depicted in figure ^. For comparison, the transition matrix 
Monte Carlo method yielded H] 

T t ~ 0.4L 28 . (17) 



We have to add here that the tunneling-time is difficult to define within the 
Wang-Landau method because it depends on the density of states which is not 
constant during the simulation. Thus r t may also depend on the interplay 
between the parameters e and f final- Thus we have chosen the typical values 
e = 0.85 and f final — 1.0 • 10~ 6 for all considered system sizes. 
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3.2 Errors at right edges 

As already mentioned in section |2.1| , one shortcoming of the single spin flip 
Wang-Landau method are the relatively large errors in the density of states at 
right edges of an energy interval over which one wants to determine the density 
of states [Q. In the N-fold way version presented here such systematic errors 
do not occur. In order to demonstrate this we have calculated the density of 
states for the first 25 levels of a L = 32 two-dimensional Ising model using single 
spin flips as well as N-fold way updates, see figure ||. This test also shows that 
the simulational effort is enormously reduced in regions where the acceptance 
rate for a spin flip is low. 

In the single spin flip Wang-Landau method systematic errors occur at right 
edges due to a boundary effect which sets in as soon as the energy interval does 
not cover all the possible energies the system may have. If the random walk 
then hits an energy level which is outside the allowed range this move is simply 
rejected and the old level is counted once more analogous to the case when no 
boundary is involved in a transition, i.e., no difference is made between the 
density of states at boundaries and those away from the boundary, whereas in 
our N-fold way algorithm the density of states at edges is treated correctly since 
it is forbidden by the definition of the flip rates (eq. (]?])) to choose a spin whose 
flipping would result in an energy outside the allowed range. 

4 Thin Ising films 

The two-dimensional Ising model is an ideal testing ground for Monte Carlo 
algorithms since it can be solved exactly and therefore permits one to verify 
simulational results directly. To show that our algorithm works also well with 
other more complicated Ising systems, we have considered the following Hamil- 
tonian 

Ti = — J SiSj — H Si — hi Si — /i2 Si with Si = ±1, (18) 
<*j> iebulk iesurfacei iesurface2 

which is essentially a three-dimensional L x L x D Ising model confined by 
two walls with fields hx,h% acting on them and periodic boundary conditions 
in the L x L planes. A sketch of this geometry is depicted in figure ||. With 
H = and hi = —hi this model exhibits the critical behavior of an interface 
localization/delocalization transition. For finite D one has a single transition 
at T C (D) which presumeably belongs to the two-dimensional Ising universality 
class. For T < T C (D) the interface is bound either to the left or the right wall 
and for T C (D) < T < T c b it is fluctuating delocalized around the center of the 
film, where T c b is the critical bulk temperature. If T satisfies T > T c b the film 
is disordered with exception of the response to the surface fields hi, hi near the 
walls. It has become clear p^-^8[ that such systems are extremely difficult to 
simulate, due to the fact that one has to deal with a much more severe slowing 
down than in simpler models, which results from the presence of a very large 
length scale £|| parallel to the surfaces in the phase with the delocalized inter- 
face which is related to the bulk correlation length via £y oc exp(Z?/4£b) pjj . 
The corresponding correlation time is then described by r oc £jf , with a dynamic 
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exponent of z m 2 [p2| . One can argue [jl6| that the asymptotic region of the two- 
dimensional Ising critical behavior is very narrow and that the true nature of the 
transition can be observed only for L>(|j. Consequently one is forced to study 
rather large linear dimensions L, depending on the thickness D of the film, since 
for D = 12, for ex amp le, one has to satisfy L 3> exp(Z)/4^6) = 26 (£h w 0.92 for 
J/k b T = 0.25 mil). In the 

case considered here (D — 6, L = 32) we have 
exp(_D/4£{,) = 5.1 and the maximum of the specific heat is still rather broad, as 
can be seen from figure 0. While for standard Ising models cluster-algorithms 
|25| , p6f have proven to be a remedy for critical slowing down their utility has to 
be tested in detail when it comes to thin Ising films with additional surface fields. 
In [^7j such a model has been studied with a ghostspin Swendsen-Wang cluster 
algorithm p8f and a significantly reduced auto-correlation time was achieved 
only for relatively small surface fields. Therefore we find it promising to test 
the applicability of our new algorithm to the model in question. For this reason 
we have performed 5 runs of a system with size 32 x 32 x 6 and competing 
surface fields hi/ J = —h 2 /J — 0.55. In order to calculate canonical averages 
of themodynamic observables we have estimated the density of states g(E) over 
the energy interval -j^ 6 [—2.83,0.2] which consists of already 185974 different 
energy levels despite the moderate size of the system. The whole interval was 
partitioned into 100 adjacent intervals with small overlaps, which are needed 
to join the density of states afterwards. For the normalization of the density 
of states we solely used the twofold degeneracy of the groundstate. The av- 
erage total number of sweeps for the complete energy range stated above was 
(1.14 ± 0.18) • 10 7 , the average total CPU time amounted to 242.4 h on a HP 
v-class, whereby a single energy interval needed (8.7 ± 0.2) ■ 10 3 s CPU time 
on average. The simulation was carried out in a parallel fashion such that the 
effective average running time was 12 h for a complete run. The simulation was 
performed with e = 0.95 and f final = (7.05 ± 0.03) • lO" 10 . 
In particular we have calculated the internal energy 

ZEg(E)e~P E 

U(T) = (E)r = E Eme ., E , (19) 

E 



the specific heat 
as well as the Gibbs free energy 



C(T) = {E2)T ~ 2 {E) - , (20) 



F(T) = -kTlnlj29( E > 

\ E 



_ y-" E \ . (21) 

V E 

and the entropy 



over a range of temperatures near Tq(D) and checked it against conventional 
Monte Carlo data obtained from a heat bath algorithm |lq] . We also considered 
canonical averages of the order parameter |m| = ^ sTf: 

E(H)e9(E)^ E 

{H)T= E E9(E)e^ > (23) 

E 
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its fourth order cumulant 



and the susceptibility 

x(T) = ^((m 2 } T -(H)|), (25) 

by sampling microcanonical averages (-Ye during the last stage of the simula- 
tions. The results are shown in figure ^ - |l2|. The value of the extrapolated 
(L — > oo) critical temperature J /ksTc{D = 6) = 0.2655 ± 0.0002 estimated in 



ref. 1 16 1 is marked by an arrow. One can see from these figures that quantities 
like the specific heat, the magnetization, the susceptibility and the fourth order 
cumulant can be obtained with an accuracy that is clearly better than that of 
the corresponding Monte Carlo data obtained with the standard single spin flip 
heat bath algorithm Jlq| , although the effort in computer time was comparable. 
Of course, for a precise analysis of the critical behaviour of the present model 
it is necessary to repeat the study for several choices of L and perform a finite 
size scaling analysis. This will be presented elsewhere |P9| . A further merit of 
the present approach is that it immediately yields precise data for the entropy 
and free energy of the model as well (figs. ||,|9|), which would be accessible from 
the standard algorithm | jl6f only by tedious techniques of "thermodynamic in- 
tegration" [|]-|| . The availability of the free energy is particularly advantageous 
in the case of first-order interface localization-delocalization transitions p8| , [l9]| , 
where the relaxation time increases exponentially with the linear dimension L 
of the system, and straightforward Monte Carlo studies would be hampered by 
hysteresis for large L (and huge statistical errors for intermediate values of L 

11). 



5 Conclusions 

In this paper, it has been shown that the efficiency of the flat histogram method 
of Wang and Landau Jl^ j can be improved significantly by combining it with 
the so-called "N-fold way" |l5| technique of sampling states of a lattice model 
via Monte Carlo methods. In particular, the problem that errors get strongly 
enhanced near the right edge of an energy interval that is sampled is elimi- 
nated, and performing systematic comparisons for the two-dimensional Ising 
model with transition matrix Monte Carlo approaches proposed by Wang and 
Swendsen |l4[ it is shown that the present algorithm is competitive in efficiency 
with techniques described in the latter study Q , but the merit of the present 
method is that it is rather straightforward to implement, and it should be imme- 
diately useful for wide classes of lattice models. Unlike cluster algorithms, there 
is no problem with the inclusion of magnetic fields. As a nontrivial example, 
we show new results of a first application to a thin Ising film with competing 
boundary fields. 



9 



Acknowledgement 



One of us (B. J. S.) received financial support from the Deutsche Forschungs- 
gemeinschaft (DFG) under grant No. Bi314/16-2. We are grateful to David P. 
Landau for numerous stimulating discussions, and for providing preprints of ref. 
|l2[ prior to publication. We are grateful to Jian-Sheng Wang and Robert H. 
Swendsen for sending a preprint of ref. pij ]. This work also received current 
support from the Bundesministerium fur Bildung und Forschung (BMBF) under 
grant No. 03N6015. 



References 

K. Binder and DW. Heermann, Monte Carlo Simulations in Statistical 
Physics. An Introduction. 3 rd ed. Springer, Berlin 1997. 

D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algo- 
rithms to Applications. Academic Press, London 1996. 

D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in 
Statistical Physics. Cambridge Univ. Press, Cambridge 2000. 

N. Metropolis, A. W. Rosenbluth, A. H. Teller and E. Teller, J. Chem. 
Phys. 21, 1032 (1953). 

G. M. Torrie and J. P. Valleau, Chem. Phys. Lett. 28, 578 (1974); G. M. 
Torrie and J. P. Valleau, J. Comp. Phys. 23, 187 (1977). 

B. A. Berg and T. Neuhaus, Phys. Lett. B 276, 249 (1991); Phys. Rev. 
Lett. 68, 9 (1992). 

A. P. Lyubartsev, A. A. Martsinovsky, S. V. Sherkunov and P. N. 
Vorontsov-Velyaminov, J. Chem. Phys. 96, 1776 (1992). 

Y. Iba, Int. J. Mod. Phys. C 12, 623 (2001). 

J. Lee, Phys. Rev. Lett. 71, 211 (1993). 

P. M. C. Oliveira, T. J. P. Pcnna and H. J. Herrmann, Eur. Phys. J. B 1, 
205 (1998). 

J.-S. Wang, Physica A 281, 147 (2000). 



F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001); fcond- 



mat/0107006 



J.-S. Wang, |cond-mat/0103318j . 



J.-S. Wang and R. H. Swendsen, |;ond-mat/0104418 . 

A. B. Bortz, M. H. Kalos and J. L. Lebowitz, J. Comput. Phys. 17, 10 
(1975). 

[16] K. Binder, D. P. Landau and A. M. Ferrenberg, Phys. Rev. Lett. 74, 298 
(1995); Phys. Rev. E 51, 2823 (1995). 



10 



[17] K. Binder, R. Evans, D. P. Landau and A. M. Ferrenbcrg, Phys. Rev. E 
53, 5023 (1996). 

[18] A. M. Ferrenberg, D. P. Landau and K. Binder, Phys. Rev. E 58, 3353 
(1998). 

[19] M. Miiller and K. Binder, Phys. Rev. E 63 021602 (2001). 

[20] P. D. Beale, Phys. Rev. Lett. 76, 78 (1996). 

[21] A. O. Perry and R. Evans, Physica A 181, 250 (1992). 

[22] P. C. Hohcnberg and B. I. Halpcrin, Rev. Mod. Phys. 49, 435 (1977). 

[23] A. J. Liu and M. E. Fisher, Physica A 156, 35 (1989). 

[24] M. Hasenbusch and K. Pinn, Physica A 192, 342 (1989). 

[25] R. H. Swendscn and J.-S. Wang, Phys. Rev. Lett. 58, 86 (1987). 

[26] U. Wolff, Phy. Rev. Lett. 62, 361 (1989). 

[27] O. Dillmann, W. Janke, M. Miiller and K. Binder, J. Chem. Phys. 114, 
5853 (2001). 

[28] J.-S. Wang, Physica A 161, 249 (1989). 

[29] B. J. Schulz, M. Miiller and K. Binder, in preparation. 



11 



20000 



15000 - 



p~ 10000 



5000 



' ^ 1 1 1 1 1 1 1 1 

10 20 30 40 50 



Figure 1: Tunneling-time rt of the two-dimensional Ising model in units of sweeps. 
Tt was fitted according to x ■ L a which yielded x ~ 0.79 and a ~ 2.58. We have used 
e = 0.85 and ff ina i ~ 1.0 ■ 10 -6 . The indicated error bars result from averages over 
30 runs (L = 8, 16, 24, 32, 36) and 10 runs (L = 40, 44, 50). 
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1000 


0.0030 


8 


(1.840 ± 0.095) • 10 Y 


0.99 


(7.0 ±0.2) • 10" lu 


1000 


0.0017 


16 


(2.19 ±0.16) ■ 10 s 


0.8 


(7.7 ±0.3) • 10"' 


100 


0.031 


16 


(9.31 ± 0.42) ■ 10 b 


0.9 


(7.0 ±0.3) ■ 10~ s 


100 


0.015 


16 


(5.7 ±0.6) ■ 10 b 


0.95 


(7.0 ±0.3) • 10~ y 


100 


0.0060 


32 


(4.2 ±0.2) ■ 10 b 


0.7 


(7.1 ±0.3) • 10"' 


100 


0.060 


32 


(8.19 ±0.46) ■ 10 b 


0.8 


(7.4 ±0.3) ■ 10~ Y 


100 


0.042 


32 


(1.065 ± 0.042) ■ 10 B 


0.85 


(7.3 ±0.3) ■ 10 -v 


100 


0.037 


32 


(7.9 ±0.6) ■ 10 B 


0.9 


(6.7 ±0.3) • 10~ M 


100 


0.016 


50 


(9.1 ±0.4) ■ 10 s 


0.7 


(7.6 ±0.2) • 10" ' 


100 


0.080 


50 


(1.371 ± 0.067) ■ 10 b 


0.85 


(7.4 ±0.2) ■ 10" ' 


100 


0.065 


50 


(9.5 ±0.8) ■ 10 B 


0.9 


(5.9 ±0.5) • 10" s 


100 


0.030 



Table 1: Simulation data for the two-dimensional Ising model. The number of sweeps, 
ffinai and the relative error s of the density of states g(E) are averages over 30 runs. 
Under checks we have listed the number of sweeps between two successive checks of 
the histogram H(E) for flatness,^ g{E) was obtained by normalizing with respect to 
the total number of states, eq. fel). e(E) for L = 5 (h old) is depicted in figure H. See 
figure H for comparison with the data provided in MM . 
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Figure 2: Relative error s(E) in the density of states g(E) for the two-dimensional 
Ising model of size L = 50. The error e(E) is an average over 30 runs. The average 
number of sweeps was (9.1 ± 0.4) • 10 5 . We have used e = 0.7 and ffi na l = (7.6 ± 
0.2) ■ 10~ 7 . The average relative error is e = 0.080. The mirror image from E/J = 
up to E/J = 5000 is not depicted. 
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Figure 3: Average error e in the density of states g(E) over the number of sweeps, as 

listed in table [j] (thick solid lines correspond to L = 50, 32, 16, 8 from top to bottom) 

in comparison with the data provided in |l4] for the same system sizes (see legend on 

right top) produced by the algorithms 

algoO: TMMC with original flat histogram flip rates. 

algo-1 exact rate: Also a TMMC with exact multi-canonical flip rates. 

two stage: Algorithm algoO followed by algo-1 which uses the estimated density of 

states from the first stage for the multi-canonical flip rates. 

algol3: Implementation of the Wang-Landau flat histogram method using single spin 
flips. 
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a i N-fold way: (4.4±0.2)-10 3 sweeps, f linsl =(7.6±0.2)-10" 
a a N-fold way: (1 .6±0.2)-10 4 sweeps, f linal =(7.4±0.2)-10" 
o o single spin flip: (1 .98+0.63) 10° sweeps, f„ nal =8.09-10 



10 



LU 



G-.-o- ...4--A— A--i - B A 



-O— -O 

-A — -A - A- — 4 A 



10"' 



-2048 -2028 -2008 -1 988 -1 968 -1 948 

E/J 



Figure 4: Relative error e(-E) in the density of states g(E) of the first 25 energy levels 
of a two-dimensional Ising model with L = 32. g(E) was obtained by normalizing with 
respect to the groundstate, e(E) is an average over 30 runs. We have used e = 0.95 
(circles and triangles) and e = 0.98 (squares). The N-fold way was approximately 35 
times faster (triangles) and 10 times faster (squares) concerning CPU time than the 
single spin flip approach. 
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Figure 5: Thin film geometry with magnetic fields hi, Y12 acting on the surfaces 
(shaded) and a field H acting on the bulk. Parallel to the surfaces periodic boundary 
conditions are imposed. 
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Figure 6: Internal energy U/N of a thin Ising film (eq. (|lij)) with L = 32 and 
D = 6. The indicated error bars resulted from an average over 5 runs. The average 
relative error of U /N in the depicted range of J/fc^Ilis 0.017%. Lines are only guides 
to the eye. The squares are data taken from ref. ]16| . The arrow marks the critical 
temperature |16| . 
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Figure 8: Gibbs free energy F/N of a thin Ising film with L = 32 and D = 6. The 
indicated error bars resulted from an average over 5 runs. The average relative error 
of C/N in the depicted range of J/k B T is 0.0011%. 



17 




Figure 9: Entropy S/N of a thin Ising film with L = 32 and D = 6. The indicated 
error bars resulted from an average over 5 runs. The average relative error of S/N in 
the depicted range of J/kgT is 0.015%. 
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Figure 10: Order parameter (|m|)T of a thin Ising film with L = 32 and D = 6. The 
indicated error bars resulted from an average over 5 runs. The average relative error 
of in the depicted range of J/kgT is 0.5%. The squares are data taken from 

ref. 
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Figure 11: Fourth order cumulant U4 of the order parameter of a thin Ising film 
with L = 32 and D = 6. The indicated error bars resulted from an average over 5 
runs. The average relative error of U4 in the depicted range of J/kgT is 6%, For 
0.26 < J/ksT < 0.27 it drops to 0.5%. The squares are data taken from ref. |16|. 
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